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SUMMARY 


The overall objective of this research program is to obtain a better 
understanding of the interaction between the solar wind and the interstellar medium 
through the use of numerical solutions of the time-dependent magnetohydrodynamic 
(MHD) equations. The simulated results have been compared with observations 
where possible and with the results from previous analytic and numerical studies. 
The primary accomplishment of this project has been the development of codes for 
2-D models in both spherical and cylindrical coordinates and the application of the 
codes to the solar wind/interstellar medium interaction. Computations have been 

carried out for both a relatively simple gas-dynamic interaction and a flow-aligned 
interstellar magnetic field. The results have been shown to compare favorably with 
models that use more approximations and to modify and extend the previous results as 
would be expected. The simulations have also been used along with a data analysis 

study to provide a quantitative estimate of the distance to the termination and bow 
shocks. Some of the specific topics that have been studied are : (1) Gasdynamic 
models of the solar wind/interstellar medium interaction, (2) Termination shock 
response to large-scale solar wind fluctuations, and (3) Distances to the termination 

shock and heliopause. The main results from each of these studies are summarized 

below. 

1. Scope of the investigation 

The interaction of the solar wind with the interstellar medium has received 

increased interest recently due to improved observations of the local interstellar 

medium and to the fact that a deep space probe may soon intersect the shock 
terminating the supersonic solar wind flow. Previous theoretical and computational 
efforts at studying the interaction have been hampered by a lack of observational 
guidance and have not been developed to the level necessary to determine the 
quantitative effects of some of the important known physical processes on the global 
interaction. When used in conjunction with present and forthcoming observations, 
the present models may assist in determining the effects of the relevant physical 

processes. In addition, they may also indirectly be able to place bounds on the 
currently unknown or speculative values of physical quantities in the local 

interstellar medium as well as the location of the termination shock and the 

heliopause. 

Time-dependent, multi-dimensional, gas dynamic and MHD codes have been 

used to study this interaction through the use of a relaxation technique in which an 

assumed non-equilibrium initial state approaches a final steady-state equilibrium. 
Although the codes incorporate numerous physical mechanisms, they can optionally 
be included in individual studies. The approach is to begin with simple cases, to 
compare earlier work and build on previous results, and to separately include effects, 
such as the interstellar and interplanetary magnetic fields, thereby providing a 
consistent basis from which to quantify the effects of a particular mechanism. 

Two separate codes for the 2-D model have been developed. In one the 

equations are solved in the r-0 plane of a cylindrical coordinate system in which the 
axial (z) axis is perpendicular to the interstellar flow direction. For the second the 
equations are solved in a spherical coordinate system. The two coordinate systems 
are used in order to allow greater flexibility for the inclusion of a magnetic field 
within the confines of a 2-D model. A variable grid spacing is incorporated in the 
radial direction to provide better resolution in the solar wind portion ( and near the 



termination shock( of the interaction. Each of the three studies carried out are 
discussed in more detail below. All of these results have been presented in 
international or national scientific meetings and published in journals. 


2. Gasdynamic models of the solar wind/interstellar medium interaction 

The first study involves the simulation of representative solar 

wind/interstellar medium interactions in both spherical and cylindrical coordinates 
when the interstellar flow is supersonic. As would be expected, both a termination 
shock and a bow shock are formed. In both coordinate systems the termination shock 
is noticeably distorted in the downstream direction. These 2-D results have been 
analyzed and compared with results from more simplified earlier studies (e.g., the gas 
dynamic study of Baranov, Space Sci. Rev., 52, 89, 1990). 

The interaction between the solar wind and the interstellar medium is 
modeled self-consistently using numerical solutions of the time-dependent 
gasdynamic equations in spherical and cylindrical coordinates. For the results 
obtained here it is assumed that the solar wind system moves through the 
surrounding medium with a supersonic velocity. After an initial (nonequilibrium) 
state has been specified, the numerical solutions follow the evolution in time until 
the interaction relaxes to a dynamic equilibrium. As would be expected, the solutions 
show the formation of a bow shock upstream of the traveling solar system to deflect 
the interstellar plasma around the cavity created by the solar wind. A termination 
shock also forms to slow and compress the solar wind plasma. For the simulation in 
spherical coordinates, the downstream portion of the termination shock reaches 
equilibrium more than three times further from the sun than the equilibrium 
distance to the termination shock on the upstream side. The results were published 
in the paper “Gasdynamic models of the solar wind/interstellar medium interaction” 
in Geophysical Research Letters, Vol. 21, 245-248, 1994. 

3. Termination shock response to large-scale solar wind fluctuations 

The second study concerns the response of the termination shock to large- 

scale fluctuations in the solar wind pressure. The analysis of data recorded by the 
Voyager 2 spacecraft indicates the presence of large-scale fluctuations in the solar 

wind ram pressure on the time scale of tens of days. The amplitude of the 
fluctuations is highly variable but often lies within a factor of 5 to 10 change from 
an average mean value of the ram pressure. Since the spacecraft has presumably 

not encountered the termination shock yet, these fluctuations should eventually 
interact with the shock and thereby play a role in determining the shock location. 

Numerical solutions of the time-dependent gas equations are used to simulate 
the response of the termination shock to fluctuations in the solar wind ram pressure 
comparable to those observed. The primary of this study is that the maximum shock 
excursion due to fluctuations is of the order of 1 AU, which is smaller than that 
predicted by other studies. Additional simulations show that the limited movement is 
due to the fact that the time scale for the termination shock response is substantially 
larger than the time scale of the fluctuations. It is also shown that the heliopause 

acts as a barrier for the fluctuations and confines them to the heliosphere. These 
results were published in the paper “Termination shock response to large-scale solar 
win fluctuations” in Journal of Geophysical Research, Vol. 99, 13,307-13,314, 1994. 



4. Study of distances to the termination shock and heliopause 


The third study is a collaboration with a data analysis effort by Dr. D. Gumett 
of the University of Iowa. We use a numerical gasdynamic simulation of an 
interplanetary shock, propagating through an equilibrium solution of the solar 
wind/interstellar medium interaction, to compute the distances to the termination 

shock and the heliopause that are consistent with the observations. 

A new heliospheric radio emission event observed by Voyagers 1 and 2 in mid- 
1992 have been interpreted by Gumett et. al. (Science, 1993] as occurring when a 

strong interplanetary shock interacts with the heliopause. The interplanetary 

shock can be related to large Forbush decreases recorded by neutron monitors at 
Earth and by cosmic ray detectors on Pioneer and Voyager spacecraft. The shock is 
thought to have originated near the Sun during a period of intense solar activity in 

late-May and early June, 1991. The time delay between the Forbush decrease at the 
Earth and the radio emission is about 408 days. The shock propagation speed as 
determined from the Pioneer and Voyager data is estimated to be within the range of 
600 and 800 km/sec, with some preference for the higher value. 

These observational results have been used as the foundation for a numerical 
study to simulate the propagation of a solar-generated shock wave through a 

dynamic equilibrium solution for the solar wind/interstellar medium interaction. A 
series of parametric studies is used to Find the equilibrium solution in which a shock 

propagating at a selected speed within the observed range (600-800 km/sec) results 
in a delay ( 408 days). The solar wind conditions at 1 AU were fixed at n = 5 cm* 3 , T = 
10 5 °K, v = 400 km/sec, and the interstellar conditions were fixed at T = 10 4o K, v=25 
km/sec. For a given interplanetary shock speed the interstellar density and 

magnetic field magnitude were varied until the above time delay was obtained. The 
termination shock is at a distance of 115 AU, and the heliopause is at 160 AU for an 

interplanetary shock speed of 800 km/sec. These distances reduce to 92 AU for the 

termination shock and 128 AU for the heliopause when the shock speed is reduced to 
600 km/sec. The results were published in the paper “Distances to the termination 
shock and heliopause from a simulation analysis of the 1992-93 heliospheric radio 
emission event” in Geophysical Research Letters, Vol. 22, 651-654, 1995. 



PUBLICATIONS 


Publications and presentations resulting from research supported by this 

NASA grant are listed below. Copies of the published papers are included in the 

appendix with this final report. 

1. R. S. Steinolfson, A numerical study of the interaction between the solar wind 
and the interstellar medium. Solar Wind VII, Goslar, Germany, 16-20, 1991. 

2. R. S. Steinolfson and T. E. Holzer, A numerical study of the interaction between 
the solar wind and the interstellar medium. XVII General Assembly of the 
European Geophysical Society, Edinburgh, Scotland, 6-10 April, 1992. 

3. R. S. Steinolfson, V. Pizzo, and T. Holzer, The interaction of the solar wind with 
the interstellar medium, EOS 73, 445, 1992. 

4. R. S. Steinolfson, and V. Pizzo, Gas-dynamic simulations of the two-shock model 
of the solar wind/interstellar medium interaction, EOS 74, 234, 1993. 

5. R. S. Steinolfson, The interaction of the solar wind with the interstellar 
medium, XXIV General Assembly of the International Union of Radio Science, 
Kyoto, Japan, 25 August - 2 September, 1993. 

6. R. S. Steinolfson, V. Pizzo, and T. Holzer, Gas dynamic models of the solar 

wind/interstellar medium interaction, Geophys. Res. Lett., 21, 245, 1994. 

7. R. S. Steinolfson, Termination shock response to large-scale solar wind 
fluctuations, J. Geophys. Res., 99, 13,307-13,314, 1994. 

8. R. S. Steinolfson and D. A. Gumett, Termination shock and heliopause distances 

estimated from coordinated simulations and data analyses, EOS, 74, 487, 1993. 

9. R. S. Steinolfson and D. A. Gumett, Distances to the termination shock and 

heliopause from a simulation analysis of the 1992-93 heliospheric radio 

emission event. Second Pioneer-Voyager Symposium, Durham, New 

Hampshire, 31 May - 3 June, 1994. 

10. R. S. Steinolfson and D. A. Gumett, Distances to the termination shock and 

heliopause from a simulation analysis of the 1992-93 heliospheric radio 

emission event, Geophys. Res., 22, 651-654, 1995. 



APPENDICES 


R. S. Steinolfson, V. Pizzo, and T. Holzer, Gas dynamic models of the solar 
wind/interstellar medium interaction, Geophys. Res. Lett., 21, 245, 1994. 

R. S. Steinolfson, Termination shock response to large-scale solar wind 
fluctuations, J. Geophys. Res., 99, 13,307-13.314, 1994. 

R. S. Steinolfson and D. A. Gumett, Distances to the termination shock and 
heliopause from a simulation analysis of the 1992-93 heliospheric radio 
emission event, Geophys. Res., 22, 651-654, 1995. 



GEOPHYSICAL RESEARCH LETTERS, VOL. 21, NO. 4, PAGES 245-248. FEBRUARY 15, 1994 


Gasdynamic models of the solar wind/interstellar medium 
interaction 

R. S. Steinolfson 

Department of Space Science, Southwest Research Institute, San Antonio, Texas 

V. J. Pizzo 1 

San Juan Capistrano Research Institute, San Juan Capistrano, California 

T. Holzer 

National Center for Atmospheric Research, High Altitude Observatory, Boulder, Colorado 


Abstract . The interaction between the solar wind and the 
interstellar medium is modeled seif-consistently using 
numerical solutions of the time-dependent gasdynamic 
equations in spherical and cylindrical coordinates. For the 
results presented here it is assumed that the solar system 
moves through the surrounding medium with a supersonic 
velocity. After an initial (nonequilibrium) state has been 
specified, the numerical solution follows the evolution in time 
until the interaction relaxes to a dynamic equilibrium. As 
would be expected the solutions show the formation of a 
bow shock upstream of the traveling solar system to deflect 
the interstellar plasma around the cavity created by the solar 
wind A termination shock also forms to slow and compress 
the solar wind plasma. For the simulation in spherical 
coordinates, the downstream portion of the termination shock 
reaches equilibrium more than three times further from the 
Sun than the equilibrium distance to the termination shock on 
the upstream side. 

Introduction 

The interaction between the solar wind and the local 
interstellar medium has been an active area of study since at 
least the early 1960s. Through a combination of theory and 
the limited observations available, researchers have been able 
to identify what arc believed to be the important physical 
processes involved in the large-scale, global interaction. The 
respective level of contribution of the various processes to the 
overall interaction remains unclear, however, due to a lack of 
both the necessary data and sufficiently developed models. 
The current state of understanding based on both theory and 
observations is summarized in several recent reviews [e.g. 
Holzer, 1989; Sucss, 1990; Baranov, 1990]. 

Although theoretical studies have not had sufficient 
observational guidance to support the development of detailed 
quantitative models consistent with the data, there is 
indication that some portion of the necessary observations 
may soon become available. A deep space probe (Pioneer or 
Voyager) may soon reach the shock that terminates the 
supersonic solar wind flow [e.g., Suess, 1990], It may also be 
PP SS ir?o t0 dctect termination shock indirectly [Kurth et 
al., 1987; Lee, 1988] through the detection of radio frequency 
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waves generated by the shock. Gumett et al. [1993] have 
interpreted some radio emission events as having been 
generated when a strong interplanetary shock interacts with 
the heliopause. Moreover, there has been an increase in the 
observational study of the interstellar medium [e.g., Frisch, 
1986; Cox and Reynolds, 1987; Davidsen, 1993; Lallemem, 
1993]. Even with these data sources, it seems clear that our 
understanding of the physical processes in the interaction and 
its complete global structure must continue to rely heavily on 
theoretical and computational models. 

We are developing computational models to study the 
3-D interaction using magnetohydrodynamic (MHD) theory. 
Our approach is to begin with the simplest models that 
contain at least some of the relevant physics and then to build 
on these models and extend them to a degree of 
sophistication not previously realized. The models are 
time-dependent* shock-capturing models so that they allow 
the capability, once a dynamic equilibrium solution has been 
computed, of simulating the temporal response of the 
equilibrium to changes in either solar wind or interstellar 
properties. We report here on equilibrium solutions obtained 
from 2-D gasdynamic models in spherical and cylindrical 
coordinates. Other gasdynamic studies that obtain features 
similar to those produced in the present simulation have been 
reported by Shima et al. [1986], Matsuda et al. [1989], and 
Baianov and Malama [1993]. However, except for the work 
by Baranov and Malama, these computations were not for the 
proper parameter regime for the solar wind/interstellar 
medium interaction. 

The interaction between the solar wind and the interstellar 
medium probably occurs within, and is largely determined by 
the properties of, the "Very Local Interstellar Medium". The 
VLISM is that part of the local interstellar medium within 
about 0.01 parsecs (2000 AU) of the Sun. As the supersonic 
solar wind flows outward toward the VLISM plasma, it 
makes a transition to subsonic flow at a termination shock 
[Parker, 1961]. This shocked solar wind interacts directly 
with the VLISM plasma at an interface (heliopause) that 
separates the two plasmas. The region of space dominated by 
the solar wind is referred to as the heliosphere. 

It is believed that the Sun moves through the VLISM at 
a speed of about 20 km/sec. Since the thermodynamic 
properties and the magnetic field strength in the VLISM arc 
not well known, it is not clear how this speed compares to 
wave propagation speeds in the VLISM, although there is 
some consensus that it is very nearly sonic. If this speed is 
supersonic, an interstellar bow shock is expected to form to 
decelerate and deflect the interstellar medium around the 
heliospheric cavity. No bow shock should form in the 
subsonic case. We will concentrate here on solutions when 
the VLISM flow speed relative to the Sun is supersonic, 
which is often referred to as the two- shock model and has 
been reviewed by Baranov [1990], 
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Solution Procedure 

The gasdynamic equations with a ratio of specific heats 
of 5/3 are solved numerically using an explicit differencing 
scheme in two spatial directions (the r, 9 directions in both 
spherical and cylindrical coordinates centered at the Sun). The 
equations and differencing scheme are discussed in more 
detail by Steinolfson [1994]. Two coordinate systems are 
considered in order to allow greater flexibility in the later 
inclusion of the magnetic field within the confines of a 2-D 
model. In cylindrical coordinates the axis is selected to be 
perpendicular to the motion of the Sun relative to the local 
interstellar medium, and all variables are assumed to be 
independent of distance parallel to the cylinder axis. In 
spherical coordinates the solution is axisymmetric about a line 
pointing into the direction of relative motion between the 
solar system and the surrounding medium. Since the velocity 
of the Sun relative to the interstellar medium (VJ is 
referenced to our Sun-centered coordinate systems, it will be 
referred to as the interstellar velocity and is assumed to be 
parallel to the 9=0° axis. 

The simulation box covers an entire hemisphere 
(0°^180°) and extends in radius from 30 to 300 AU. Three 
of the physical variables (radial flow speed v r , electron 
density n e , and pressure p) are symmetric across the poles 
while the fourth is antisymmetric (meridional flow speed v 0 ). 
The supersonic solar wind is specified at 1 AU, and an 
adiabatic atmosphere with constant radial flow speed is used 
to compute values at the inner radial boundary for the 
computation at 30 AU. Interstellar plasma conditions and 
velocity are specified at the outer radial boundary for 
O°£0£9O°, and zero-order extrapolation along the local flow 
direction is used to obtain values at the outer radial boundary 
for 9O O <0£18O°. For fixed physical quantities in the solar 
wind and interstellar medium, the numerical computation 
continues until a dynamical equilibrium is obtained. 


Numerical Results 

Spherical Coordinates . The physical quantities in the solar 
wind at 1 AU for this simulation are v r =300 km/sec, v e =0, 
n e =3 cm' , and temperature T=10 3 *K. The conditions at 1 AU 
are used to normalize all quantities. The flow speed is 




referenced to the sound speed (52.45 km/sec), and the 
thermodynamic quantities are referenced to their values at 1 
AU. The inflow interstellar conditions are V il =24.88 km/sec, 
n =0.2 cm \ and T=l(r °K (M^l.5). 

The radial grid spacing is fixed at 1 AU, and the angular 
grid spacing is fixed at 1°, for a grid of 271x181. The initial 
(nonaquilibrium) state for the simulation is constructed by 
specifying adiabatic solar wind conditions within 30-40 AU, 
a termination shock at 40 AU, interstellar conditions beyond 
80 AU, and a linear fit between the shocked solar wind and 
interstellar plasmas. The solution continues for 44,000 time 
cycles, which represents a total physical evolution time of 
1.81x10 s days or 500 years. Based on the sound speed in the 
interstellar medium (16.59 km/scc), this is the time for a 
disturbance to cross the 600 AU numerical box more than 
twice. 

The velocity streamlines and contour plots of the 
thermodynamic quantities in the dynamic equilibrium state 
are shown in Figure 1. Symmetries, of course, could be used 
to reflect the solution about the lower boundary. The 
approximate locations of the bow and termination shocks 
have been superimposed on the computer- generated plot of 
the velocity streamlines. In terms of the commonly used 
nomenclature for shocks in the solar wind, the bow shock is 
a forward shock, and the termination shock is a reverse 
shock. The heliopause is easily identified in the figure as the 
division between interstellar and solar streamlines and has 
been extended to the pole at 0° by the dashed line. 

The spatial structure for several of the physical quantities 
from the inner boundary out to 200 AU along the radial line 
at 0° is presented in Figure 2. Although the shocks and 
heliopause are spread over several grid points (as usual for 
simulations such as this), the quantities have the expected 
variation across the discontinuities. The changes in quantities 
across the shocks agree to within a few percent of the values 
computed from the Rankine-Hugoniot equations. 

Although not shown here, the entropy has been verified 
to be constant in the regions between the discontinuities in 
Figure 2. The use of an adiabatic solar wind has the effect of 
allowing the Mach number to become very large (due to the 
monotonically decreasing temperature) at large heliocentric 
distances. Just upstream of the termination shock in Figure 2, 
the Mach number has a value of about 68, while it is about 
58 at the inner computational boundary and 5.7 at 1 AU. 




Fig. 1. Thermodynamic quantities and velocity streamlines in the relaxed dynamic equilibrium. 
The interstellar flow (relative to the solar system) is from left to right. The shock locations drawn 
on the streamline plot were determined from the contour plots of the thermodynamic variables. 
For the thermodynamic quantities in the contour plots, the value used to construct the plots is 
(Q-QJ/Qi„ where Q* represents the interstellar value. A contour level of zero represents the 
interstellar value, and values greater (smaller) than the interstellar value are represented by solid 
(dashed) contours. The density contours range from -0.9 to 1.1 in increments of 0.1, pressure 
from -0.8 to 2.6 in increments of 0.2, and temperature from 8 to 96 in increments of 8. 


Stcinoifson ct al.: Gasdynamic Models 


247 



portion of the termination shock in this case is ejected out of 
the numerical box. Preliminary studies indicate that the 
equilibrium distance of the termination shock from the Sun 
near the 180 pole increases as the interstellar Mach number 
increases. The interstellar Mach number for this run is simply 
so large (M*=12) that the equilibrium distance is beyond the 
outer computational boundary. Other than the different 
geometry, this solution is qualitatively similar to the previous 
one. The intention is to use this 2-D cylindrical model as a 
test-bed for evaluating the effects of magnetic field 
orientations that would not be possible with the above 2-D 
model in spherical geometry. Although the absolute values of 
the magnetic effects may differ from those in spherical 
geometry, the relative effects should be similar and should 
provide useful insight into the interpretation of results from 
eventual 3-D simulations. 

Discussion 


Fig. 2. The variation of physical values along the radial line 
at 0^ in the equilibrium solution. The pressure has been 
multiplied by 10 and the temperature by lu 2 for the indicated 
scale. Although the "discontinuities" are spread over several 
grid points, the change in quantities across them agrees with 
analytic expressions. 


As might be anticipated from the results in Figure 1 
(pamcularily the streamlines), a large shear flow develops 
across the heliopause away from the pole at 0°. This is 
illustrated in Figure 3, in which the variations of the 
velocities along the line A-B marked on the plot of the 
steamlines in Figure 1 are shown. The velocities are the 
components perpendicular (1) and parallel (|) to line A-B 
defining the normal to the heliopause. The shaded region 
indicates the approximate numerical width of the heliopause. 
Despite this large shear, there is no indication of the 
development of a Kelvin-Hclmholtz instability, although this 
instability has b<ren interpreted as forming at the heliopause 
in other simulations of the interaction of a stellar wind with 
the interstellar medium (Matsuda, 1989]. 

Except for the region near the pole at 18 (f, the solar wind 
plasma must turn through a large angle in order to flow out 
downstream along the heliopause, as indicated by the 
streamlines in Figure 1. Some part of the turning occurs at 
the termination shock through a reduction of the flow speed 
parallel to the shock normal. Note that the flat portion of the 
termination shock near the 180° pole turns the flow away 
from the pole. Since the termination shock in this model is a 
gasdynamic shock, the component of the speed perpendicular 
to the shock normal must be conserved across the shock. This 
is no longer required for MHD shocks, so some important 
differences would be expected when the magnetic field is 
included. MHD simulations by Washimi [1993] indicate that 
the magnetic field may produce other significant differences 
from the gasdynamic results presented here. 

Cylindrical Coordinates. The physical quantities in the solar 
wind at L f° r casc m v r =60 km/sec, n c =0.25 cm' 3 , 
and T=l(r °K. As for the above example, the thermodynamic 
quantities are referenced to their values at 1 AU, and the flow 
speed is referenced to the sound speed at 1 AU. The 
interstellar values are V u =200 km/sec, ry=8xl0* 5 cm' 3 , and 
T=*10r °K (Mj,=12). The selected physical values are clearly 
not physically realistic and have been chosen to obtain a 
two-shock solution within a few hundred AU of the Sun. Due 
to the cylindrical geometry, the solar wind density falls off 
with distance as r r rather than as r 2 as it would for spherical 
geometry. The inner boundary is at 10 AU, and the outer 
boundary is at 450 AU. The grid spacings are constant at 3 
AU and 3° for a grid of 150x60. 

The pressure is shown in a gray-scale plot in Figure 4 
along with representative streamlines. The downstream 


Two examples of numerical simulations in spherical and 
cylindrical coordinate systems for the interaction of the solar 
wind with the interstellar medium (for the case when the 
solar system is traveling supersonically relative to the 
surrounding medium) show at least one qualitatively similar 
result. That is, the distance of the termination shock from the 
Sun on the downstream side (in the solar system wake) of the 
moving solar system is substantially larger than the 
corresponding distance on the upstream side. For the example 
in spherical coordinates, the distance from the Sun to the 
termination shock in the direction opposite to that of the 
relative motion between the solar system and the interstellar 
medium is more than three times larger than the Sun to 
termination shock distance in the direction of relative motion. 
The downstream portion of the termination shock moves out 
of the computational box for the example in cylindrical 
coordinates. It should be noted that in simulations for which 
the solar system moves subsonically with respect to the 
interstellar medium for the same parametric values as those 
used here (except for the interstellar wind speed) the distance 
to the termination shock is only slightly larger downstream 
than upstream [Steinolfson, 1994]. 

The solar wind plasma obviously cannot continue to 
expand outward adiabatically on the downstream side since 
the thermal pressure would become far too small to maintain 
lateral force equilibrium with the surrounding interstellar 
plasma. We suggest that the situation is qualitatively similar 
to that for supersonic flow in a diverging nozzle exiting into 
a region with a larger back pressure than the thermal pressure 
in the flowing plasma [e.g., Shapiro, 1953]. The analogy is 



0 20 40 60 SO 

DISTANCE (AU) 


Fig. 3. The variation of the flow speeds with distance along 
the line A-B indicated on Figure 1 where the distance is 
referenced to point A on the line. The flow speeds are the 
values perpendicular and parallel to the line. 
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Fig- 4. Streamlines and gray-scale plots of the thermal 
pressure for the dynamic equilibrium solution in cylindrical 
coordinates. 


only approximate, of course, since the solution on the 
downstream side must maintain an equilibrium (and must be 
self-consistent) with the solution on the upstream side. In the 
case of the supersonic nozzle, the equilibrium position of the 
shock moves downstream as the back pressure decreases, 
when the interstellar velocity is subsonic, the back pressure 
idenncal to the interstellar pressure [Steinolfson, 
," 4 1- However, as shown in Figure 1, the back pressure in 
the present case in considerably less than the interstellar 
pressure, and, as a result, the termination shock comes to 
equilibrium turther downstream than in the subsonic study. 
Additional simulations over a wide range of parameter space 
analogy^ C ° esta ^ sfl the generality of this suggested 

TTie solar wind/interstellar medium interaction computed 
in the present work and shown in Figure 1 agrees 
qualitatively with the results obtained by Shima et al. [1986- 
shown in their Figure 2], although they were considering a 
more general stellar example and did not use physical values 
in the proper parametric regime for the solar case. Shima et 
al. used a similar shock-capturing numerical scheme, but their 
numerical resolution was substantially coarser than for the 
simulations presented here. Matsuda et al. [1989] repeated the 
Shima et ai. simulation utilizing much finer numerical 
resolution and found that the Kclvin-Helmholtz instability 
developed on the heliopause (see Figure 4 in Matsuda et al.). 

f k in “T. r™, 2 ? result not indicate the formation of 
the Kelvin- Helmholtz instability, Matsuda et al. naturally 
attributed its appearance in their simulation to the better 
numerical resolution. The numerical resolution in the present 
simulations is even finer than that in the Matsuda et al. study 
and yet there is no indication of the formation of any 
instability on the hcliopausc (see Figures 1 and 4). It is not 
clear whether this is due to the different parametric regimes 
tor the two computations, or if it is a result of the numerics 
in one of the simulations. It does seem curious that the 
instability in the Matsuda et al. study develops near the pole 
directed mto the interstellar wind (0=0°) , where the velocity 
shearis very small, and is then convccted downstream. 

The present results also agree with those obtained in a 
recent study by Baranov and Malama [1993], which is for 
parametric values applicable to the solar case although they 
are somewhat different than those used here. Baranov and 
Malama used a discontinuity-fitting numerical scheme in 
which the steady-state equations are solved between the 
whose locations are adjusted until an 
equilibrium is reached. The formation of a time -dependent 
instability, such as the Kelvin-Helmholtz instability, could not 
be studied with this numerical scheme. Consequently, the 
results of the present time-dependent simulations support the 


use of a discontinuity-fitting scheme, at least for the solar 
case. The complexity (both numerical and physical) of the 
interaction of a stellar wind with the interstellar medium 
dictates that various numerical approaches along with 
analytical work are necessary in order to establish the 
physical validity of the model interaction. 
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Termination shock response to large-scale solar wind fluctuations 
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Abstract. The analysis of data recorded by the Voyager 2 spacecraft indicates the presence 
of large-scale fluctuations in the solar wind ram pressure on the time scale of tens of days. 
The amplitude of the fluctuations is highly variable but often lies within a factor of 5 to 10 
change from an average or mean value of the ram pressure. Since the spacecraft has 
presumably not encountered the termination shock yet, these fluctuations should eventually 
interact with the shock and thereby play a role in determining the shock location. Numerical 
solutions of the time-dependent gasdynamic equations are used to simulate the response of the 
termination shock to fluctuations in the solar wind ram pressure comparable to those observed. 
The primary result of this study is that the maximum shock excursion due to the fluctuations 
is of the order of 1 AU, which is much smaller than that predicted by other studies. 
Additional simulations show that the limited movement is due to the fact that the time scale 
for the termination shock response is substantially larger than the time scale of the 
fluctuations. It is also shown that the heliopause acts as a barrier for the fluctuations and 
confines them to the heliosphere. 


1. Introduction 

As the solar wind plasma flows outward from the Sun, it 
eventually interacts with the surrounding interstellar plasma, 
and the two plasmas relax to a dynamic pressure equilibrium. 
The formation of a cavity in the interstellar medium by the 
solar wind plasma was initially suggested by Davis [1955]. 
For a pressure equilibrium between the solar and interstellar 
plasmas to be reached, the solar wind makes a transition 
from supersonic to subsonic flow at a termination shock 
[ Parker , 1961]. The termination shock probably occurs 
within the range of 50 to 100 astronomical units (AU) from 
the Sun [e.g., Suess, 1990], At some distance beyond the 
termination shock, the solar plasma interacts directly with the 
interstellar plasma at an interface separating the two plasmas. 
This interlace is referred to as the heliopause, and the region 
containing solar plasma is referred to as the heliosphere. 

It is generally believed that the solar system travels 
through the interstellar medium at a speed of approximately 
20 km s‘ . The thermodynamic conditions of the interstellar 
medium are not well known, but if the speed of the solar 
system through this medium is supersonic, a bow shock 
forms beyond the heliopause to slow and deflect the interstel- 
lar plasma. Baranov [1990] has reviewed this two-shock 
model of the solar wind/interstellar medium interaction. A 
bow shock would not be expected to form when the speed of 
the solar system relative to that of the surrounding medium 
is subsonic. 

As summarized in recent reviews on this topic [e.g.. 
Holier, 1989; Suess, 1990], the genera] structure of the 
interaction of the solar wind with the interstellar medium in 
terms ot the bulk or large-scale properties is fairly well 
understood in a qualitative sense tor steady plasma flows. 

Copyright 1994 by the American Geophysical Union. 
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There are, of course, several physical processes of a more 
detailed or localized nature occurring in the two plasmas that 
are not well understood themselves (e.g., interstellar neutral 
gas, charge exchange, cosmic rays). Their effect on the 
overall interaction is also not well understood. 

The analysis of data from the plasma science experiment 
on Voyager 2 as the spacecraft traveled from 1 to 40.4 AU 
indicates that there are large-scale fluctuations in the solar 
wind ram pressure occurring on the time scale of tens of 
days [Belcher et al., 1993]. Data from Pioneers 10 and 11 
also show the presence of large variations in the ram 
pressure [ Barnes , 1990]. Since the ram pressure plays a 
major role in controlling the location of the termination 
shock, it seems logical that large-scale fluctuations in the 
ram pressure may produce relatively large excursions of the 
termination shock. Belcher et al. [1993] used a kinematic 
model in which they limited the termination shock speed to 
200 km s' 1 and estimated that the location may vary by as 
much as 15 AU in response to the observed ram pressure 
fluctuations. This excursion amplitude is similar to that 
computed in other one-dimensional analyses [Barnes, 1993; 
Suess , 1993]. If there are large-scale movements of the 
termination shock in response to large-scale fluctuations in 
the solar wind ram pressure, then a spacecraft that reaches 
the vicinity of the shock may cross it more than once. 

The objective of the present study is to use numerical 
solutions of the gasdynamic equations to examine the effect 
that large-scale fluctuations with temporal and spauai scales 
similar to those observed have on the average or mean 
behavior of the interaction between the solar wind and 
interstellar plasmas and particularly on the location of the 
termination shock. The study is performed for the case of 
subsonic relative motion between the solar system and the 
interstellar medium (for reasons discussed in section 3.1). A 
dynamic equilibrium solution of the interaction is obtained 
first, and fluctuations in the solar wind ram pressure are then 
generated to perturb the equilibrium. 
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2. Equations and Model 

The gasdynamic equations used in this study can be 
written in the following nondimensional form; 

v V- (pv) = 0, da) 

C7t 



where d/dt is the total derivative, and the physical variables 
P» v, and p are the density, plasma velocity, and thermal 
pressure, respectively. The density p =n/m pt where n e is the 
electron number density and m p is the proton mass. All 
variables have been made dimensionless by normalizing the 
thermodynamic quantities to the initial values at 1 AU and 
the velocity to the initial sound speed at 1 AU (a Q ). Dis- 
tance is referenced to 1 AU (r^, and time is normalized by 
l o~ r J a o' Th e only parameter that enters directly into the 
above equations is the specific heat ratio y, which we take to 
be 5/3. Note that the potentially important effects of 
collisions with neutral species, resonance charge exchange, 
and cosmic ray pressure have not been included. Donohue 
and Zonk [1993] used a gasdynamic model of a localized 
section of the termination shock to show that acceleration of 
the anomalous component and/or galactic cosmic rays may 
significantly influence the nature, structure, and motion of 
the termination shock. It remains to include the physics of 
their study into a global model such as that used here. In 
addition, since this is a gasdynamic model, the magnetic 
fields of the interstellar medium and of the solar wind are 
not included. The neglect of the magneuc field in the solar 
wind should have a minimal effect, since the plasma beta 
(ratio of thermal pressure to magnetic pressure) is generally 
large throughout the solar wind. It has been suggested that 
the magnetic field may provide the largest contribution to the 
total pressure in the interstellar medium [e.g., Suess y 1990], 
but the physical nature of the interstellar pressure should not 
significantly affect the present results. 

Equations (la) - (lc) are solved in the r-0 plane of a 
spherical coordinate system aligned such that the pole of the 
coordinate system at 9=0° points into the direction of motion 
of the solar system relative to that of the interstellar medium. 
The solution is independent of the azimuthal angle § and 
thus is axisymmetric about the poles. The coordinate system 
is centered at the Sun (r=0) and remains fixed with respect 
to the Sun. Consequently, the relative motion between the 
solar system and the interstellar medium is represented by a 
velocity of the interstellar medium (V ]S ) in our coordinate 
system. The simulation box extends from 0° to 180° in 0 
and from 10 AU to 400 AU. The radial grid spacing is 
constant at 1.4 AU, and the angular grid spacing is constant 
at 1°, which gives a grid of 280x18 1. 


The equations are solved numerically using a two-step 
Lax-Wendroff differencing scheme similar to that given by 
Richtmyer and Morton [1967], with second-order accuracy. 
As is usual for such schemes, an additional smoothing must 
be performed on the values from the Lax-Wendroff step in 
order to remove high-frequency oscillations and the over- 
shoots and undershoots created near ideal discontinuities such 
as shocks. A smoothing term suggested by Lapidus [1967] 
was used, as was a flux-corrected-transpon term [Book et al., 
1975]. The two smoothers produced qualitatively the same 
results when they both worked. The Lapidus term, though, 
was considerably more robust in handling the nonuniform- 
ities and the substantial range in physical variables in the 
present problem and was used for the results presented here. 
The flux-corrected-transport correction required more adjust- 
ment of the numerical constants and for some parametric 
regimes produced questionable results. 

An explicit, shock-capturing scheme, as used here, is not 
necessary for computing the initial dynamic equilibrium 
solution into which the fluctuations are introduced, although 
it is essential in order to simulate the proper temporal 
response to the fluctuations. The initial equilibrium could 
equally well have been computed using a shock-fitting 
scheme with either an explicit or an implicit method used to 
solve for the solution between the discontinuities [e.g., 
Baranov and Malama , 1993]. It might be expected that the 
addition of a semi-implicit term in the momentum equation 
would permit the use of a larger time step during the 
relaxation. However, even with such a term, the time step 
is still restricted so that plasma cannot be convected farther 
than a grid spacing in a single time step. With the small 
grid spacing and the large flow speeds in the solar wind, the 
time step for a semi-implicit method would be only about 
12% larger that for the explicit method used here. Such a 
small increase in the time step certainly would not compen- 
sate for the additional CPU time necessary to solve the 
momentum equation implicitly in a semi-implicit method. 

Three of the physical variables (radial flow speed v r 
density, and pressure) are symmetric at the angular bound- 
aries, while the fourth is antisymmetric (meridional flow 
speed v e ). The supersonic solar wind is specified at 1 AU, 
and an adiabatic atmosphere with constant radial flow speed 
is used to compute values at the inner computational bound- 
ary at 10 AU. Since this is a supersonic-inflow boundary, all 
physical quantities can be specified on it. The interstellar 
density and velocity are held fixed at the outer radial 
boundary for 0°<9<90 o , and zero-order extrapolation along 
the local flow direction is used to update the pressure. Since 
this portion of the outer boundary is a subsonic-inflow 
boundary, not all of the physical quantities can be specified 
there. Zero-order extrapolation along the local flow direction 
is used to obtain all values at the outer radial (subsonic- 
outflow) boundary for 90°<9<180 o 

The simulation is performed in two separate steps. A 
dynamic equilibrium is obtained first, and then a fluctuating 
solar wind flow speed is imposed on this equilibrium at the 
inner radial boundary. The simulation to compute the 
dynamic equilibrium is initiated by specifying adiabatic solar 
wind conditions within 10-40 AU, a termination shock at 40 
AU, interstellar conditions beyond 80 AU, and a linear fit 
between the shocked solar wind and the interstellar plasmas. 
Starting with this nonequilibrium state, the time-dependent 
simulation continues until the solution relaxes to a dynamic 
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equilibrium. The soiar wind fluctuations are then generated 
in the equilibrium solution by temporal oscillations in the 
solar wind speed at the inner radial boundary. The thermo- 
dynamic conditions at the inner boundary remain unchanged 
during the flow speed oscillations. 

3. Numerical Results 

3.1. Dynamic Equilibrium Solution 

The spatial variation of the thermodynamic quantities and 
the streamlines after the solution has relaxed to a dynamic 
equilibrium are shown in Figure 1 . The physical quantities 
in the solar wind at 1 AU for this simulation are v r =300 km 

s' 1 , V e =0, rt e =3 cm' 3 , and temperature ^lO^K. The refer- 
ence sound speed is then 52.45 km s**, and the reference or 
characteristic time (g is 33.2 days. The fixed, inflow 
interstellar conditions are V u =13.28 km s* 1 and n e =0.2 cm' 3 . 
The inflow interstellar temperature is initially set to 10 4o K 
(A/ i3 =0.8), but this temperature will change slightly as the 
pressure adjusts at the inflow boundary. The relaxation 
required 120,000 time cycles, which represents a total 
physical evolution time of 4.94xl0 5 days, or 1353 years. 
Based on the sound speed in the interstellar medium (16.59 
km s or 3.48 AU/yr), this is time for a disturbance to cross 
the 800 AU numerical box almost six times. This total 
evolved time should not be regarded as physically relevant; 
it is simply an indicator of the time required for the solution 
to adjust to the specified start-up conditions. However, 


simulated changes in this equilibrium solution due to changes 
in the boundary conditions should provide realistic estimates 
of the actual response time of the solar wind/interstellar 
medium interaction. 

The interstellar flow (in our Sun-centered coordinate 
system) is from left to right in Figure 1. The locations of 
the termination shock and the heliopause are indicated on the 
streamline plot. The rapid changes in variables at the 
termination shock appear as clustered contours on the 
pressure and temperature plots. The shock does not appear 
on the density plot because of the linear scaling used for the 
contours and the fact that the density at the shock is less than 
the minimum contour. There is a density jump by about a 
factor of 4, however, since the termination shock is a strong 
shock. For this example with subsonic interstellar flow, the 
termination shock is nearly spherical and does not have the 
characteristic bullet shape typical of solutions for a superson- 
ic interstellar flow [e.g., Baranov and Malama , 1993; 
Steinolfson et a/., 1994]. The formation of a more spherical 
termination shock is one reason a subsonic interstellar flow 
solution is used for the initial state. However, since the 
present study concentrates on examining the fluctuations 
along the 0° pole, the shape of the downstream portion of the 
termination shock should have virtually no effect on the 
results. 

3.2. Response to Large-Scale Fluctuations 

A sinusoidal temporal variation in the solar wind radial 
flow speed is now specified at the inner radial boundary in 



Figure 1 . Thermodynamic quantities and velocity streamlines in the relaxed dynamic equilibrium solution. 
The radial range extends from the inner computation boundary at 10 AU to 400 AU. The interstellar flow 
(relative to the soiar system) is from left to right. The termination shock location drawn on the streamline 
plot was determined from the contour plots of the thermodynamic variables. For the thermodynamic 
quantities, the value used to construct the plots is {Q-QQIQ^ where (9 1S represents the interstellar value. 
A contour level of zero represents the interstellar value, and values greater (smaller) than the interstellar 
value are represented by solid (dashed) contours. The density contours range from -0.96 to 0.24 in 
increments ot 0.08, pressure ranges from -0.8 to 0.5 in increments of 0.1, and temperature ranees from 10 
to 100 in increments of 10. 
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order to generate fluctuations in the soiar wind ram pressure 
(pv 2 ). The speed oscillation is' maintained throughout the 
numerical computation. Although the change in the boundary 
conditions occurs only in the flow speed (the thermodynamic 
quantities are held Fixed at the inner boundary), fluctuations 
in all the physical quantities are generated in the soiar wind. 
The same variation is applied at all angular locations, so the 
physical situation that is modeled is one in which the flow 
speed increases and decreases uniformly at all locations on 
the solar surface. The amplitude of the flow speed oscilla- 
tion is 100 km s' 1 , so the speed varies between 200 and 400 
km s' 1 at the boundary. The fluctuation period is taken to be 
180 days. This is somewhat larger than the observed large- 
scale variations in ram pressure occurring on time scales of 
tens of days reported by Belcher et at. [1993], but it will be 
sufficient to demonstrate the essential physics involved in the 
response of the interaction to the fluctuations. 

The fluctuations in ram pressure at the inner boundary ( 10 
AU) and at two other spatial locations upstream of the 
termination shock are shown in Figure 2. The ram pressure 
fluctuations at the inner boundary (10 AU) are due to 
changes in just the flow speed, but those at the other two 
locations involve changes in both the flow speed and density. 
The fluctuations are shown for only the first 2500 days of 



TIME (DAYS) 

Figure 2. Temporal fluctuations in the solar wind ram 
pressure at three locations along the pole at 0=0° upstream 
of the termination shock. The fluctuations are shown for 
only the first 2500 days but remain unchanged throughout 
the remainder of the simulation. 



Figure 3. Temporal fluctuations in the locations of the 
termination shock along the poles at 0=0° and 0=180° and of 
the heliopause along the pole at 0=0°. The fluctuations are 
shown for only the first 2500 days but remain unchanged 
throughout the remainder of the simulation. 

the simulation, but they remain unchanged for the total 
computation time of 7800 days. The termination shock in 
the equilibrium solution is located at about 80 AU; hence the 
fluctuations at 70 AU are only 10 AU upstream of the shock. 
The average value of the simulated ram pressure is less than 
the observed average value in Figure 3 of Belcher et al. 
[1993] owing to the lower- than -observed solar wind speed 
used in the simulations. This lower speed was used so that 
the termination shock would be located fairly close to the 
Sun when reasonable (anticipated) values were used for 
interstellar conditions. The magnitudes of the simulated ram 
pressure fluctuations in Figure 2 are at the low end of the 
observed magnitudes in Figure 3 of Belcher ei al. [1993] and 
are similar to those of the fluctuations during solar minimum 
periods. 

The temporal responses of the termmauon shock at the 
poles (0=0°, 180°) and of the heliopause at the 0° pole are 
shown in Figure 3. The period of the response is the same 
as that of the applied speed oscillation (180 days). The 
amplitude of the response represents a total movement of the 
termination shock of about 1 AU. It should be noted that the 
amplitude of the response shown here is virtually identical to 
that when the applied flow speed oscillation amplitude at the 
inner boundary is either reduced to 75 km s 1 or increased to 
125 km s' 1 . This indicates that the amplitude of the termi- 
nation shock response may have saturated and that it would 
not increase significantly with further increases in the 
amplitude of the ram pressure fluctuations. If the flow speed 
oscillation amplitude is decreased to 50 km/s, there is a 
noticable reduction (about 20%) in the oscillation amplitude. 
The oscillation amplitude at which the termination shock 
response saturates naturally depends on the oscillation period, 
but this dependence is not investigated here. 

Although the simulated termination shock does not move 
as much as might be anticipated in response to the ram 
pressure fluctuations, the amplitudes of the spatial oscilla- 
tions are significant, as shown in Figure 4. This figure 
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Figure 4. Spatial fluctuations in the density and flow speed along the poles (at 8=0°. 180°) at the end of 
the simulation (dashed curves) superimposed on the dynamic equilibrium state (solid curves). 


compares the spatial structures of the density and radial flow 
speed along the poles in the dynamic equilibrium solution 
( solid curves) with the structures alter the speed oscillation 
has been applied for 7800 days (dashed curves). It is evident 
from this figure that even though only a tlow speed oscilla- 
tion was applied at the inner boundary, large oscillations in 
both tlow speed and density are generated in the solar wind. 
Notice trom the plots along the 0° pole that the fluctuations 
do not get through the heliopause into the interstellar 
medium. This is more apparent in Figure 5, in which the 
contour levels have been tailored to bring out fluctuations in 
the thermal pressure beyond the termination shock. The 
heliopause effectively acts as a barrier to contine the oscilla- 
tions to the solar wind plasma. 

The simulated response amplitude is substantially less than 
the 10 AU (or larger) movement of the termination shock 
due to ram pressure fluctuations predicted by more approxi- 
mate theories [Barnes, 1993; Suess, 19931. It is also less 
'ban the shock movement estimated in the study by Belcher 


et ai [1993], in which the termination shock location was 
computed using observed ram pressure fluctuations with a 
200 km s l limit on the shock speed. There are at least two 
reasons why the results of the studies cited above for the 
response amplitude of the termination shock differ from the 
present results. 

First, the earlier studies are all one-dimensional and 
effectively assume that an isolated section of the termination 
shock responds independently of the rest of the termination 
shock. This, of course, does not happen in practice, and the 
entire termination shock responds collectively to solar wind 
ram pressure tluctuations. In the extreme case, one can 
envision situations in which one section of the termination 
shock is being subjected to a maximum in a ram pressure 
fluctuation while an adjoining section, not too far removed 
in angular distance from the first but at a different radial 
distance trom the Sun, must respond to a minimum in the 
fluctuation. The first section will tend to move outward, and 
'he second one will tend to move inward. Although this 
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Figure 5. Contour plots of the thermal pressure illustrating confinement of the fluctuations to the 
heliosphere. The heavy dashed curve marks the location of the heliopause. The fluctuations within the 
termination shock do not appear in this figure, since the pressure within the termination shock is below the 
minimum contour level. 


etfect must be accommodated to a certain extent for all 
fluctuations, it becomes more important as the period of the 
fluctuations decreases. A second, much more important 
reason for the difference between the present results and 
those in earlier studies probably has to do with the more 
rapid response times either calculated or assumed in the 
earlier studies compared to those computed in the present 
study. This is now discussed in more detail. 

3.3. Termination shock response time 

As an estimate of the time it takes the termination shock 
to respond to changes in the solar wind quantities, we rapidly 
increase the solar wind flow speed at the inner boundary 


from 300 km s' 1 in the equilibrium solution to 400 km s 1 , 
hold the speed constant at 400 km s 1 at the inner boundary! 
and follow the evolution in time until the interaction ap- 
proaches a new dynamic equilibrium. The flow speed is 
increased over 45 time cycles, which represent about 90 
days, or half the period of the fluctuations used in the above 
study. The temporal response of the termination shock 
locations at the poles and of the heliopause at the 0° pole are 
given in Figure 6. The termination shock and heliopause 
gradually approach their new equilibrium position. 

The termination shock at the 0° pole initially moves 
outward rapidly (between about 300 and 700 days in Figure 
6) when it is first impacted by the solar wind disturbance. 
The shock at this location then rebounds (between about 700 



TIME (DAYS) 

Figure 6. Temporal responses of the termination shock locations at the poles and of the heliopause at the 
pole : as the solar wind speed in the dynamic equilibrium solution in Figure 1 is increased over 90 days 
tram 300 km s' 1 to 400 km s 1 and then at that level. 
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and 1000 days) and then monotonicaliy approaches an 
equilibrium position. The average outward speed of the 
termination shock along the 0° pole during the initial rapid 
movement, approximately 73 km s' 1 , is followed by an 
average outward speed from 1000 to 6000 days of 7.2 km 
s . Both of these speeds are smaller than the maximum 
shock speed computed or estimated in earlier studies [Bel- 
cher et al , 1993; Barnes , 1993; Suess , 1993], During the 
90-day half-period of the fluctuation used above, the termina- 
tion shock would move 3.8 AU if traveling at 73 km s' 1 and 
0.38 AU if traveling at 7.2 km s' 1 . These distances are 
smaller than the excursions predicted in the earlier studies 
and are, respectively, larger and smaller than the 1 AU 
movement reported on above (section 3.2). The average 
speed of the termination shock during the computed 1 AU 
fluctuations is about 20 km s' 1 (1 AU in 90 days). This 
study also shows that if a single increase in ram pressure of 
infinite duration is input, as Barnes [1993] did, the termina- 
tion shock would then have time to reach larger speeds. 

The termination shock at the 180° pole also experiences 
an initial overshoot, but the heliopause does not. The motion 
ol the termination shock along the 180° pole following the 
initial overshoot in Figure 6 illustrates an important feature 
ot the relaxation simulation. Once the downstream (from 
90° to 180°) portion of the termination shock reaches the 
vicinity of its equilibrium position, it oscillates over a very 
long period. This particular computation was not continued 
beyond the time shown in Figure 6, so it is not known 
whether the oscillation amplitude decreases with time. The 
simulation to determine the initial equilibrium solution (sec- 
tion 3.1) was continued long enough to determine that the 
oscillation amplitude was slowly decreasing with time. The 
spatial variation of the physical quantities and the streamlines 
for a 400 km s' 1 solar wind speed are not shown here. 
However, other than the larger distances to the termination 
shock and the heliopause, there are no qualitative differences 
between that solution and the solution shown in Figure 1 for 
a 300 km s 1 solar wind speed. 


4. Discussion 

The response of the termination shock to large-scale 
fluctuations in the solar wind ram pressure is studied using 
numerical solutions of the gasdynamic equations. The first 
step in this study is to compute a dynamic equilibrium 
solution tor the interaction between the solar wind and the 
interstellar medium by using a numerical relaxation proce- 
dure. The ram pressure fluctuations are then generated in the 
solar wind by temporal flow speed oscillations at the inner 
radial boundary for the computation. The subsequent solar 
wind ram pressure fluctuations are comparable in amplitude 
(near the observed minimum) and time scale (near the 
observed maximum) to those observed by Voyager 2 [Belch- 
er et al.. 1993]. 

The maximum movement of the termination shock in 
response to the ram pressure fluctuations is about 1 AU for 
a 100 km s 1 flow speed oscillation at the inner boundary. 
Although only an oscillation in the flow speed is used as the 
driving mechanism, the solar wind disturbance consists of 
fluctuations in all physical quantities, including the density. 
When the driving flow speed oscillation is reduced to 75 km 
s , the termination shock response is almost identical to that 


for the iarger-magmtude speed oscillations reported on here 
(which, in fact, is the same as that for 125 km s' 1 oscilla- 
tions), indicating that the response has saturated and would 
remain the same for even larger amplitude fluctuations. As 
a result, the fact that the magnitude of the ram pressure 
fluctuations are near the low end of those observed should 
not make any significant difference. That is, the termination 
shock response would be the same as that computed here 
even if the oscillation amplitude were increased to the 
maximum observed value. A reduction to 50 km s' 1 , 
however, reduces the termination shock movement some- 
what. It should also be noted that the time scale for the 
fluctuations used here (180 days) is somewhat larger than 
observed. Consequently, the predicted 1-AU excursion 
should be considered an upper limit which would most likely 
be reduced for shorter-period fluctuations. One additional 
caution that should be mentioned is that the ambient or back- 
ground solar wind flow speed of 300 km s' 1 used in the 
study is lower than typical observed values. The effect of a 
larger ambient speed should be investigated, though one 
would anticipate that the deviation from the ambient value 
would have more of a role in determining the amplitude of 
the termination shock fluctuation. 

The magnitude of the termination shock excursion 
computed here (1 AU) is substantially less than the 10 to 15 
AU excursions predicted by other methods [Belcher et al., 
1993; Bames y 1993; Suess , 1993]. These other methods are 
more approximate in the sense that they are one-dimensional 
and thus consider only the movement of an isolated section 
of the termination shock rather than the collective response 
of the entire termination shock. The previous methods also 
do not consider the complete interaction region between the 
solar wind and the interstellar medium. One possible 
explanation for this difference is that the numerical differenc- 
ing scheme used here contains inherent dissipation, as does 
any such numerical scheme, which may act to reduce the 
termination shock response time. This effect has been 
minimized in the present study by not adding any additional 
dissipation. Along the same lines, it should be emphasized 
that the differencing scheme is explicit, so all relevant time 
scales are included in the evolution. 

The response of the termination shock may also be 
restricted by the fact that it is spread over several grid points 
in the simulation. This effect is difficult to quantify, but it 
should be noted that the spatial scale of the oscillations is 
considerably larger than the numerical width of the termina- 
tion shock, which should tend to minimize the effect. A 
much more likely explanation for at least the major portion 
of the difference is that the termination shock simply does 
not respond as rapidly as either assumed or calculated in the 
earlier studies. The response time was simulated in the 
present study by computing the time needed for the interac- 
tion to reach a new equilibrium following a fixed change in 
the solar wind properties and was shown to be much larger 
than that used in the previous studies. However, given the 
numerical approximations and uncertainties, the actual 
termination shock response will probably be bounded by the 
values predicted here and those from the earlier work. 
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Abstract. A new heliospheric radio emission event 
observed by Voyagers 1 and 2 in mid-1992 is believed 
to have been produced by the interaction of an inter- 
planetary shock with the heliopause. The shock is 
thought to have originated near the Sun during a period 
of intense solar activity in late-May and early-June, 

1991. The observed travel time of the shock to the 
heliopause is 408 days; the initial speed is estimated to 
be between 600 and 800 km/s. We use a numerical 
gasdynamic simulation of an interplanetary shock, 
propagating through an equilibrium solution of the solar 
wind/interstellar medium inter action, to compute the 
distances to the termination shock and the heliopause 
that are consistent with these observations. For a shock 
speed of 600 km/s, the termination shock is located at 
92 AU, and the heliopause is located at 128 AU, These 
distances increase to 112 AU and 156 AU when the 
shock speed is increased to 800 km/s. 

Introduction 

Pioneers 10 and 11 and Voyagers 1 and 2 are on their 
way to explore the interaction of the solar wind with the 
interstellar medium. Pioneers 10 and 11 are currently at 
heliocentric distances of 58.6 and 39.8 AU (astronomical 
units), and Voyagers 1 and 2 are at 54.5 and 41.9 AU (all 
distances are as of January 1, 1994). None have yet crossed 
the termination shock, where the solar wind is expected to 
become subsonic, or the heliopause, the boundary between 
the solar wind and the interstellar plasma. Because of the 
exploratory nature of these missions, there is considerable 
interest in obtaining improved estimates of the distances to 
these boundaries. Most previous estimates have been ob- 
tained by assuming a total pressure for the interstellar 
medium and computing boundary locations from pressure 
balance considerations. Unfortunately, relatively little is 
known about the various contributions to the total pressure 
of the interstellar medium in the vicinity of the Sun. Based 
on current best estimates, the distance to the termination 
shock is believed to be between 70 and 100 AU, and the 
distance to the heliopause between 100 and 150 AU [e.g.: 
Holzer, 1989; Baranov, 1990; Suess, 1990; Suess and 
Nemey, 1994]. A strong, new heliospheric radio emission 
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event was recently detected by the plasma wave instruments 
on the Voyagers 1 and 2 spacecraft [Gumett et al., 1993]. In 
their analysis of this event, Gumett et al. proposed that the 
radio emission is produced by the interaction of an interplan- 
etary shock with the heliopause. Making rather simple 
assumptions about the speed of the shock, Gumett et al. 
estimated that the heliocentric radial distance to the helio- 
pause is from 117 to 177 AU. Here we use a numerical 
simulation to more accurately model the propagation of the 
interplanetary shock and thereby obtain an improved estimate 
of the distance to the heliopause (and termination shock). To 
obtain a solution that agrees with the observations, the total 
pressure of the interstellar plasma is varied until the travel 
time of the interplanetary shock through the corresponding 
equilibrium solution for the solar wind/interstellar medium 
interaction agrees with the observed travel time. With this 
approach, the interstellar pressure is a result of our study 
rather than an assumed quantity. These numerical results 
depend, of course, on the validity of the interplanetary shock 
model proposed by Gumett et al. [1993]. 

The 1992-93 Heliospheric Radio Emission Event 

Beginning in early July 1992, the plasma wave instru- 
ments on the Voyagers 1 and 2 spacecraft (located at 49.0 
AU and 37.6 AU, respectively) detected strong, new helio- 
spheric radio emission around 2 to 3 kHz [Gumett et al., 
1993]. This event reached peak intensity in early December 

1992, and declined to near the receiver noise level by mid- 

1993. At peak intensity the total power radiated was estimat- 
ed to be at least 10 13 Watts. Gumett et al. [1993] proposed 
that the radio emission was produced by the interaction of an 
interplanetary shock with the heliopause. The interplanetary 
shock is believed to have originated at the Sun during a 
period of intense solar activity in late-May and early-June, 
1991. The coronal mass ejections and associated interplane- 
tary disturbances associated with the iate-May/early-June 
solar events are believed to have merged in the outer 
heliosphere into a single, quasi-sphericai interplanetary shock 
that was subsequently detected by various instruments on 
Pioneers 10 and 11 and Voyagers 1 and 2. The propagation 
speed of the interplanetary shock has been variously estimat- 
ed to be 820 ± 45 km/s [Van Allen and Fillius, 1992; Van 
Allen, 1993], 600 to 800 km/s [Webber and Lockwood, 
1993], and 550 km/s [Belcher et al., 1993]. The uncertainties 
are due in part to the difficulty of identifying the exact onset 
times at the Sun, and in part to the spread in the arrival 
times at the various spacecraft. Since the shocks may 
decelerate, their speeds at large distances from the Sun may 
be overestimated. 
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The travel time of the shock from the onset of the period 
ot intense solar activity on day 145, 1991, to the arrival at 
the heliopause on day 188, 1992, is 408 days, or 1.12 years 
[Gumett et al, 1993], The travel time is also subject to 
various uncertainties, and in fact, McDonaid et al. [19931 
have suggested that the interaction may have been triggered 
by an earlier period of solar activity in March 1991 . Because 
of the uncertainties in the propagation speed and travel time 
of the interplanetary shock, for the purpose of this paper the 
original values reported in Gumett et al. [1993] will be used; 
i.e„ 600 to 800 km/s and 408 days. 

Numerical Procedure 

The gasdynamic model used here is the same as that used 
by Steinolfson [1994], For the present study two separate 
simulations are required. In the first simulation, a dynamic 
equilibrium solution is obtained by numerically computing 
the relaxation from an initial nonequilibrium state with 
specified solar wind conditions at 1 AU and specified 
conditions at the interstellar inflow boundary. In the second 
simulation, an interplanetary shock with a speed of either 
600 or 800 km/s is introduced into the equilibrium solution 
near the Sun and allowed to propagate outward through the 
solar wind and into the interstellar medium. 

For the equilibrium solution, the physical quantities in the 
solar wind at 1 AU are held fixed at the following values- 

? e f U L Cm ' J ’ tem P erature T e =10 3o K. radial velocity 
(v r ) E -400 km/s, and theta velocity (v 9 ) E =0. The solar wind 
Mach number at 1 AU is M^.6. The quantities in the 
interstellar plasma that are held fixed are the temperature 
U.-10 K and the interstellar flow speed V, =24.88 km/s 
These values give a supersonic interstellar flow with a Mach 
number of Mi,= 1.5. A trial value is selected for the interstel- 
ar density n,,, and a dynamic equilibrium solution is comput- 
ed. The final interstellar density is determined iteratively by 
requiring that an interplanetary shock with a speed of either 

600 or 800 km/s propagate from the Sun to the heliopause in 
408 days. 

Sun-fixed spherical coordinates are used in the simula- 
tions. The simulation box extends from 30 AU to 550 AU 
and from 0“ to 180° in 0, with the 0=0“ pole directed into the 
interstellar flow. The solution is assumed to be axisymmetric 
about the poles so there is no variation with the azimuthal 
angle <() and the computation becomes two-dimensional A 
coordinate transformation allows the radial grid spacing to 
vary monotonically from a minimum of 0.85 AU at the inner 
boundary to a maximum of 3.15 AU at the outer boundary 
ihe angular grid spacing is constant at 1°, which gives a grid 
ot 335x181. The equations are solved numerically using a 
second-order explicit scheme with a high-frequency filter 

Three ot the physical variables (v,, n, and pressure p) are 
symmetric at the poles while the fourth is antisymmetric (v 0 ) 
An adiabatic solar wind with constant flow speed is used to 
obtain boundary values at the inner radial boundary at 30 AU 
tram the specified solar wind conditions at 1 AU Since this 
is a supersonic-inflow boundary, ail physical quantities can 
be specified on it. The outer radial boundary for O°<0<9O° is 
also a supersonic-inflow boundary. Zero-order extrapolation 
along the local flow direction is used to obtain values at the 
outer radial [outflow] boundary tor 9O°<0<18O O 

After a dynamic equilibrium solution has been computed, 
an interplanetary shock propagating at a speed of either 600 
or 800 km/s in separate simulations is generated at the inner 
radial boundary. The shock is produced by changing the 


boundary conditions at 30 AU to the values given by the 
shock jump conditions tor a shock at the inner boundary 
traveling at the selected speed. These revised boundary 
conditions are then maintained for the duration of the 
simulation. 

Numerical Results 

The dynamic equilibrium solution in which a shock 
initiated at the inner boundary with a speed of 600 km/s 
takes approximately 408 days to travel from the Sun to the 
heliopause is shown in Figure 1. A constant 600 km/s shock 
speed from the Sun to the inner computational boundary at 
30 AU is assumed. The pressure is referenced to the inter- 
stellar value. Due to the linear spacing used for the contours, 
the spatial variation of the pressure within the termination 
shock is not shown in this representation. Although not 
shown here, the entropy has been verified to be constant 
along the streamlines between the discontinuities. The 
pressure and other thermodynamic plots have been used as 
a guide to draw the approximate locations of the shocks and 
heliopause on the computer-generated streamline plot. The 
termination shock is elongated in the downstream direction 
and has the bullet shape characteristic of the equilibrium 
solution for a supersonic interstellar flow [e.g., Baranov and 
Malama, 1993; Steinolfson et al., 1994]. 

The pressure increases from the bow shock to the 
heliopause near the 0=0° pole to form a high pressure region 
just ahead of the heliopause. Since the flow is adiabatic 
between the bow shock and the heliopause, a density pile-up 




Figure 1 . Thermal pressure contours and velocity stream- 
lines in the dynamic equilibrium solution for a solar wind 
shock speed of 600 km/s. The interstellar flow (relative to 
the solar system) is from left to right. The shock locations 
plotted on the streamline plot were determined from the 
contour plots of the thermodynamic variables. The value 
used to construct the plots for the pressure is (Q-QJ/Q ia 
where Q lt represents the interstellar value. Values greater 
(smaller) than the interstellar value are indicated by solid 
(dashed) contours. The pressure contours range from -0.8 to 
2.4 in increments of 0.2. 


region is also formed ahead of the heliopause. Such a density 
pile-up has been used by Gumett et al. [19931 to explain the 
increase in the radio emission frequency with increasing time 
as the shock propagates through the region beyond the 
heliopause. 

The spatial variation of the thermodynamic quantities and 
the flow speed along the radial line at the 0° pole are given 
in Figure 2. The thermodynamic quantities are normalized to 
their values at 1 AU, and the flow speed is referenced to the 
sound speed at 1 AU in this representation. The shaded 
regions identify the approximate locations of the shocks and 
the heliopause. Although the shocks and heliopause are 
spread over several grid points, the changes in physical 
quantities across them are consistent with the jump condi- 
tions across shocks and contact surfaces. 

The temporal variations of the radial location (dashed 
curves) and velocity (solid curves) along the 0° pole of an 
interplanetary shock propagating through the above equilibri- 
um solution are shown in Figure 3. The velocity is positive 
when the shock motion is away from the Sun. In this 
example the initial shock speed (in the laboratory frame) at 
the inner boundary is 600 km/s, slowing to 540 kirn's when 
the shock reaches and interacts with the termination shock. 
This interaction results in the generation of an inner shock 
which is the new (perturbed) termination shock, and an outer 
shock which is the continuation of the interplanetary shock. 
The perturbed termination shock initially moves outward at 
about 190 km/s. It overshoots its final equilibrium position, 
and after a strongly damped oscillation comes to rest at a 
distance slightly greater than 107 AU. 

During the interaction with the termination shock, the 
interplanetary shock quickly decelerates from 540 km/s to 
450 km/s and then slows by another 40 km/s before reaching 
the heliopause. At the heliopause the shock again decelerates 
quickly from 410 km/s, eventually approaching an asymptot- 
ic speed of about 40 km/s as it propagates into the interstel- 
lar medium. During the time period represented by the light 
lines, the numerical logic used to track the interplanetary 
shock does not work well due to the complex shock-shock 
interactions. 

As can be seen in Figure 3, the interaction of the inter- 
planetary shock with the heliopause occurs over several days. 



RADIAL DISTANCE (AU) 

Figure 2. The variation of physical quantities along the 
radial line at 0=0° in the equilibrium solution. The pressure 
has been multiplied by 10 and the temperature by 10* 3 for the 
indicated scale. The thermodynamic quantities are referenced 
to their values at 1 AU using Q/Q E , where Q E is the value at 
1 AU. The flow speed is divided by the sound speed at 1 
AU. Negauve values ot the flow speed are directed toward 
the Sun. 



Figure 3, The temporal variation of the velocity (solid 
curves) and location (dashed curves) along the 0° pole of the 
interplanetary shock prior to and after its encounter with the 
termination shock and of the termination shock after it has 
been set in motion by the interaction with the interplanetary 
shock. The horizontal dashed line indicates a zero shock 
velocity. 


To compare with the observations, a time toward the end of 
the interaction is taken as the onset time of the observed 
radio emission. To obtain the total time delay from the 
injection of the interplanetary shock near the Sun to the onset 
of the radio emission, the time necessary for the shock to 
travel from the Sun to the inner boundary at 30 AU must be 
added to the time given in Figure 3. For a shock speed of 
600 km/s, it takes about 65 days for the shock to traverse 30 
AU. Consequently, the total travel time corresponds to the 
observed value of 408 days. As seen in Figure 2, the 
termination shock is centered at about 92 AU and the 
heliopause is at 128 AU along the 0° pole for the dynamic 
equilibrium solution. The interstellar plasma density has a 
value of ^,=0.14 cm* 3 for this solution. 

A study identical to that described above was also 
performed using an interplanetary shock speed of 800 km/s. 
In this case the termination shock is located at 1 1 2 AU and 
the heliopause at 156 AU. Since the interplanetary shock is 
now traveling faster, the termination shock and the helio- 
pause must be located farther from the Sun to give the 
observed 408-day propagation time. The interstellar density 
required to produce this new equilibrium configuration is 
correspondingly lower, ^,=0.09 cm 3 . Since there is no 
qualitative difference from the previous 600 km/s solution, 
illustrations analogous to those shown above are not given 
for the 800 km/s solution. 

The interstellar plasma densities obtained for the above 
two equilibrium solutions are probably unphysically large 
due to the omission of the interstellar magnetic field pres- 
sure. This omission is not expected to adversely affect the 
results of the study performed here, especially if the field 
orientation is such that tension forces on the heliopause are 
minimal (such as for a How-aligned magnetic field). Let us 
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tor now assume that only the total interstellar pressure is 
important in the pressure balance, without regard for how the 
pressure is obtained. Then, it we further assume that the 
interstellar magnetic field provides some of the computed 
total pressure and that the interstellar density has a more 
physically realistic value, the magnetic field strength neces- 
sary to provide this pressure contribution can be easily 
computed. For our values of T 1( and V IS , the sum of the 
thermal and dynamic pressures becomes UxlO' 11 dynes/- 
cm . The magnetic field strength in nT necessary to maintain 
this total pressure with an interstellar density (in cm' J ) less 
than that used in the simulations can be written as 

B - ■ . 

where (nj,^ is the value from the simulation and (n if ) ob# is 
the value from observations. 

For a shock speed of 800 km/s, the density from the 
simulation is (nJ^O.09 cm' 3 . If the density in the interstel- 
lar medium is (nj ob- =0.04 cm' 3 , as suggested by Gumett et 
al. [1993], then B^O.41 nT, whereas if (nJ^O.Ol cm 3 . 
B it =0.5 1 nT. The corresponding values for the 600 km/s 
shock speed, for which CnJ tinJ =0. 14 cm 3 , are B, =0.57 nT for 
^Job,=0.04 cm J and B 1 =0.65 nT for (nj obg =0.01 cm' 3 . The 
difference between the minimum [0.41 nT] and the maxi- 
mum [0.65 nT] values for the magnetic field strength is quite 
small. 

Discussion 

Gumett et al. [1993] have proposed that the 1992-93 
heliospheric radio emission event is caused by the interaction 
of an interplanetary shock with the heliopause. The numeri- 
cal simulations presented in this paper demonstrate that for 
an interplanetary shock speed of 600 km/s the observed 
travel time of 408 days can be obtained if the interstellar 
total pressure is such that the termination shock in the 
equilibrium solution is located at 92 AU and the heliopause 
is located at 128 AU. If the shock speed is increased to 800 
km/s, these distances increase to 112 AU and 156 AU. 

If the distances computed in the present numerical study 
are taken to be realistic estimates, the earliest any spacecraft 
would encounter the termination shock would be in 2004 by 
Voyager 1 [based on the spacecraft trajectories in Figure 2(b) 
from Suess, 1990], Voyager 2 would reach the termination 
shock in 2008 and Pioneer 1 1 in 2014. The above encounter 
times are based on the estimates using the 600 km/s shock 
speed. If the results for the 800 km/s shock speed are used, 
the encounter with the termination shock increases to 2009 
for Voyager 1, 2014 for Voyager 2, and 2018 for Pioneer 11. 
These encounter times are underestimates since the distances 
to the termination shock found in this study are those along 
the direction of relative motion between the Sun and the 
interstellar medium, which, as seen in Figure 1, are the 
minimum distances. All the spacecraft are traveling at a 
relatively large angle to this direction and, as a result, would 
require more time to reach the termination shock. 

One factor that might compromise the estimated distances 
to the termination shock and the heliopause determined from 
this study is the fact that the interplanetary shock is initiated 
at 30 AU rather than at the surface of the Sun so any shock 
deceleration within 30 AU is not modeled self-consistently. 
In fact, the shock is assumed to have a constant speed from 
the Sun to 30 AU. This difficulty is mitigated somewhat in 
that the time the shock spends in traversing 30 AU is about 


15 % of the total time required for the shock to reach the 
heliopause. It is certainly conceivable that our initial inter- 
planetary shock speeds at 30 AU may be too large due to 
shock deceleration. A lower interplanetary shock speed 
would reduce the distances to the termination shock and the 
heiiopause. The spread of the shocks and heliopause over a 
few grids also makes it hard to select a time when the 
interactions between the interplanetary shock and the 
termination shock and heliopause actually occur. In addition, 
the interplanetary magnetic field strength, not included in our 
gasdynamic study, increases from the termination shock to 
the heiiopause. This increases the fast magnetosomc speed 
and thereby modifies the interplanetary shock speed. Despite 
these qualifications and others of less importance, the 
distances estimated here should be accurate to within 
something less than 10 AU. 
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